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Abstract. 

We study the collective states of interacting non-Abelian anyons that emerge in 
Kitaev's honeycomb lattice model. Vortex-vortex interactions are shown to lead to the 
lifting of the topological degeneracy and the energy is discovered to exhibit oscillations 
that are consistent with Majorana fermions being localized at vortex cores. We 
show how to construct states corresponding to the fusion channel degrees of freedom 
and obtain the energy gaps characterizing the stability of the topological low energy 
spectrum. To study the collective behavior of many vortices, we introduce an effective 
lattice model of Majorana fermions. We find necessary conditions for it to approximate 
the spectrum of the honeycomb lattice model and show that bi-partite interactions are 
responsible for the degeneracy lifting also in many vortex systems. 
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1. Introduction 

One of the main trends in contemporary condensed matter physics is the study of 
topological order. Apart from their fundamental interest, the motivation derives from 
the discovery of topological quantum computation [U |2] . The quasiparticle excitations 
of topologically ordered systems are anyons that can be employed to perform fault 
resilient quantum information processing [31 It is of interest to find out what are the 
sufficient conditions for topological order to exist, which systems are useful for quantum 
computation and how anyons emerge. While the research concentrated initially on 
the experimentally challenging fractional quantum Hall states or abstract spin lattice 
constructions, the recently discovered topological insulators have given rise to novel 
constructions of anyon supporting systems |5]. The experimental accessability of the 
latter has led to the theory and detection of Majorana fermions becoming the subject 
of intense research. 

Majorana fermions are the simplest possible non-Abelian anyons [6]. While they 
are not universal for quantum computation [71 [8] , their realization and the subsequent 
demonstration of topological information processing is considered a significant stepping 
stone towards a full scale implementation. Physical systems that are predicted to 
give rise to Majorana fermions include p-wave superconductors [9], Kitaev's honeycomb 
lattice model [10] and its generalizations [HI |12l [131 El US] , the Moore- Read state in 
fractional quantum Hall systems [I6] and lately also various superconductor/insulator 
interfaces [171 [TBI . Experiments have been performed on fractional quantum Hall 
states [20] as well as on topological insulator interfaces [5], but the existence of Majorana 
fermions remains still ambiguous. For the detection of Majorana fermions, and for their 
future applications, it is crucial to understand how their anyonic properties are manifest 
in the microscopic physics. The essential property to all topological quantum computing 
schemes is the existence of non-local topological degrees of freedom. For Majorana 
fermions they appear as fusion degrees of freedom that describe the different locally 
indistinguishable ways the anyons can behave collectively. However, these states can 
become locally distinguishable due to tunnehng processes that lead to the degeneracy 
being lifted [2Tj. This has been studied in Kitaev's honeycomb lattice model [22], 
in the Moore- Read state [23] as well as in p-wave superconductors [Ml ES, [26]. It 
has also recently been discovered that in many anyon systems the tunneling processes 
can not only lift the degeneracy, but can even drive transitions between topological 
phases [271 ISHl [29] . The physical conditions under which this occurs, however, depend 
on the microscopies of a specific model. Only when the tunneling is exponentially 
suppressed, the topological phase is stable and its anyonic excitations can be employed 
for topologically protected quantum computing. 

In the first part of the paper we analyze the degeneracy lifting in the non- 
Abelian phase of Kitaev's honeycomb lattice model [lOj. By simulating continuous 
vortex transport [30], we extend the vortex- vortex interaction analysis of [22j. We 
find the anticipated oscillations in the energy splitting show how they depend 
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on the microscopies of the model and characterize the stabihty of the topological 
low energy spectrum by obtaining the relevant energy gaps. In the second part we 
introduce an effective Majorana fermion lattice model, an extension of the tight-binding 
model considered in [31], to study the degeneracy lifting in many vortex systems. 
Simultaneous interactions between vortices lead in general to a hybridization of the 
states corresponding to fusion degrees of freedom [32]. We find necessary conditions, i.e. 
relevant effective flux sectors, for the Majorana model to approximate this subspectrum. 
The flux sectors are interpreted to arise from the vortices of the underlying model, 
which provides a direct connection between the microscopic and the effective models. 
Furthermore, we find that the effective tunneling amplitudes are well approximated by 
the energy splitting due to vortex-vortex interactions. This confirms that the picture of 
freely tunneling Majorana fermions applies also in many vortex systems. 

This paper is organized as follows. In Section II we review the solution of the 
honeycomb lattice model using Majorana fermionization. We show how it gives rise to 
an equivalence between coupling and vortex configurations and how this can be employed 
to simulate continuous vortex transport. In Section III we analyze the oscillating 
degeneracy lifting of the fusion degrees of freedom when the separation between two 
vortices is varied. We obtain the physical parameter dependence of both the oscillations 
and the energy gaps characterizing the topological low energy spectrum. In Section IV 
we introduce a lattice model of Majorana fermions and employ it to study the collective 
behavior of up to eight vortices. 

2. The honeycomb lattice model 

In this section we briefly review Kitaev's honeycomb lattice model and its solution by 
mapping it to free Majorana fermions. For a more rigorous treatment we refer to the 
original work [10] and to the subsequent developments [22[ |33l |3ll [35l |36l EH |38] . We 
show how to relate the manipulation of the vortices to the manipulation of the model's 
physical parameters. This enables the simulation of continuous vortex transport |30j . 
which is later used to study the physics of the anyonic vortices. 

2.1. Solution by mapping to free Majorana fermions 

Kitaev's model [lOj, comprises of spin-1/2 particles residing on the vertices of a 
honeycomb lattice. The spins interact according to the Hamiltonian 

a=x,y,z {i,j)^a— links {i,j,k) 

where J^j are real nearest neighbour couplings on links of type a specified by their 
orientation (see Figure[T](a)). The second term is an effective magnetic field of magnitude 
K, which explicitly breaks time-reversal invariance. The sum in this term runs over all 
next to nearest neighbour triplets as described in [22] . 
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Figure 1. The conventions for tlie lioneycomb lattice model, (a) The links are oriented 
as shown and labeled x, y oi z depending on their orientation, (b) A single plaquette p 
with its six sites enumerated, (c) The oriented six next to nearest neighbour couplings 
per plaquette originating from the three spin term in the Hamiltonian |10j . (d) The 
normalized lattice basis vectors and rij^. 

The model can be mapped to a free fermion problem when the spins are represented 
by Majorana fermions pi]]. This procedure reduces the Hamiltonian to the quadratic 
form 

if = - ^ hijCiCj, hij = 2JijUij + 2K^ UikUjk- (2) 

i,j k 

Here c, = cj are Majorana fermions that satisfy {cj, Cj} = 26ij. The first term describes 
their hopping between nearest neighbour sites while the second term describes next to 
nearest neighbour hopping between sites shown in Figure [Tt^c). The Hermitian operators 
Uij = ulj live on the links of the honeycomb lattice and should be understood as a Z2 
gauge field. Since [uij,H] = 0, they have no dynamics. The gauge field interpretation 
follows from the local constraint | \E') = | \1') that the eigenstates | of the original 
Hamiltonian ([T]) have to satisfy for all sites i . The operator Di anti-commutes with 
the three link operators acting on vertex i and thus it can be interpreted as a local gauge 
transformation. The Hamiltonian is gauge invariant. 

At the heart of the exact solvability of the model are the local symmetries 
[wp,H] = 0, where, as illustrated in Figure [U^b), Wp = 0"iC"2'^3C"4C5C"6 Hermitian 
plaquette operators associated with every plaquette p. In the fermionized picture these 
become 

Wp = M12M32M34M54M56M16, [wp, Di] = 0, (3) 

that in the gauge theory language can be understood as gauge invariant Wilson loop 
operators. Their eigenvalues Wp = —1 are interpreted as having a vr-fiux vortex on 
plaquette p. The pattern of all plaquette operator eigenvalues w = {wp} is referred 
to as a particular vortex sector, that is created by fixing the gauge, i.e. choosing the 
pattern of the link operator eigenvalues u = {uij}. The vortex sectors label the physical 
sectors of the model. 

We consider finite systems of 2LxLy spins with and Ly plaquettes in directions n^; 
and My (see Figure [H^d)). Due to particle-hole symmetry in the problem, diagonalization 
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of ([2]) in vortex sector w yields in general the double spectrum 

i=l 

where hi = Yl'j=i'^[^T]j'^j} ^'^^ some complex coefficients [Q:f]j, are complex fermionic 
modes. When the diagonalization is performed numerically, ±Ef and cxf follow from 
solving the eigenvalue problem ih^ct^ = —Efctf. Here h'^ is the real kernel matrix ([2]) 
when restricted to vortex sector w. 

2.2. Gauge/coupling configuration equivalence and vortex transport 

Our aim is to study the physics of the vortices by considering the spectral evolution as a 
function of the vortex sectors w. While the sectors form a discrete set, we can effectively 
interpolate between them, i.e. simulate the transport the vortices, as follows. As can be 
seen from ([2]), the local value Uij of the gauge field appears always uniquely paired with 
a local coupling J^-. This implies that while the vortex sector depends on the gauge 
through the plaquette operators, ([3]), from the point of view of the Hamiltonian, (|2]), we 
can regard u as the sign configuration of the couplings. It follows that the Hamiltonians 
H'^ and H^' on two different vortex sectors can always be connected by flipping the 
signs of some of the couplings, i.e. one can always find some coupling configurations J 
and J' such that H'^iJ) = H^' {J') % Therefore, by manipulating the local couplings we 
can interpolate the spectrum continuously between different vortex sectors and thereby 
study the spectral evolution when vortices are created, transported and annihilated. 
When this is performed adiabatically, the vortex sector w does not strictly speaking 
change, but in our simulation the spectrum will evolve as if it had. This same method 
has been employed to evaluate also the non-Abelian statistics of the vortices [30]. 

From now on we adopt a perspective where the system has been initially prepared 
in the ground state belonging to the vortex- free sector [39]. The vortex sectors w are 
viewed as coupling configurations J = {Jij} of equal amplitude, but with locally varying 
signs. To transport a vortex between two plaquettes or to create/annihilate a pair of 
vortices on them, one needs then to change the sign of the physical coupling Jj^ on the 
shared link. This can be performed continuously by tuning Jij — — t- —Jij which 
we simulate by changing the sign of the coupling in S steps such that at step s its 
value is given by Jijil — As 5* becomes large, this approximates well a continuous 
process. The resulting transport protocol is illustrated in Figures [2l^a)-(d). Treating the 
vortex sectors w and the coupling configurations J on equal footing is also experimentally 
motivated. Given sufficient site addressability, the local control of the couplings through 
external laser fields is how one could perform vortex creation and adiabatic transport 
in the proposed optical lattice experiments [lOl E] . 

X Formally one should also manipulate the local values of K. However, when K is small compared to 
the couplings Jij , the magnitude of the three spin term does not significantly affect the physics of the 
vortices and thus controlling the couplings Jij is sufficient. 

§ Introducing a complex phase is not possible, because the couplings have to be real at all times. 
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Figure 2. A protocol for simulating continuous vortex transport, (a) Initially the 
coupling configuration is chosen such that = — I on the links crossed by the dashed 
line, while Jij = I on all other links. This corresponds to a vortex on the left plaquette. 
(b) Consider changing the coupling on the link in the middle from J^j = 1 to Jij = — 1 
in S steps of size ^ . At step s its value is J-j = 1 — ^ , which we interpret as the vortex 
occupying a location away from the plaquette center, (c) When Jfj = 0, the Wilson 
loop operator is defined only on the composite plaquette, and the vortex is interpreted 
being located in the middle, (d) Finally, as Jfj — — 1, the vortex moves smoothly to 
the plaquette on the right. 

2. 3. System of interest 

To study the physics of the anyonic vortices, we employ a large finite toroidal system 
of 1200 spins defined by Lx = 30 and Ly = 20 and consider the system initially in the 
vortex-free sector by setting Uij = 1 on all links (any other gauge giving rise to Wp = 1 
on all plaquettes would also do). While this sector of the honeycomb lattice model 
supports several topological phases where the vortices can behave as either Abelian or 
non-Abelian anyons [lOj, we concentrate here only on the latter. It is characterized by 
Chern number z/ = ±1, which implies that the vortices obey Ising anyon statistics. This 
phase emerges when K ^ and the couplings violate the inequahties Ja> J13 + J7 for all 
all permutations of a, /3, 7 = x, y, z. Unless otherwise noted, we achieve this by setting 
J = = = 1 on all links. This amounts to considering the system in the middle 
of the non-Abelian phase in the sense that the spectral gap is maximized as a function 
of the couplings Ja- For J < 1 the system is closer to the phase boundary located at 
J = |. On top of this background we then imprint different vortex configurations by 
locally changing the couplings as discussed above. 

We retain the magnitude K of the three spin term as a free parameter. Its physical 
range depends on the way it is introduced in the proposed optical lattice experiments 
[40l |4T| . If it is introduced by applying a weak Zeeman field [10], only magnitudes of 
up to i^' ~ 10~^ can be tolerated before the topological phase breaks down [12]. On the 
other hand, by engineering it directly one can in principle achieve larger values [13]. 

3. Vortex-vortex interactions and the fusion degrees of freedom 

In this section we show that in the presence of vortices the spectrum exhibits modes 
whose energy oscillates and converges to zero with the vortex separation. These modes 
are argued to correspond to the fusion degrees of freedom of the non-Abelian anyonic 
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Figure 3. (a) Two vortex configuration parametrized by the vortex separation dg 
(picture not on scale for the — 30, Ly = 20 system). The blue squares denote the 
location of the vortices, which are created by setting = — 1 on the links crossed by 
the dashed blue line while = 1 on all other links, (b) The particle-hole symmetric 
mode spectrum Q as a function of ds- The two dashed red lines are the oscillating 
fusion modes with energies ±5''= , while the solid black lines are the vortex separation 
independent free fermion modes with energy gap A/. The plot is for K = 0.05. 

vortices. We find the microscopic dependence of tlie oscillations and show that they 
are consistent with Majorana fermions being localized at the vortex cores. Finally, we 
outline the topological low energy spectrum. 

3.1. Oscillating vortex-vortex interactions 

When two vortices are separated linearly as illustrated in Figure ^a), the spanned 
configurations are conveniently parametrized by a continuous vortex separation defined 
as ds = d + ^ . Here d denotes the number of links of equal coupling strength between 
the vortices and 1 < s < S is the step when the instantaneous coupling strength on 
the {d + l)th link is 1 — ^. At integer values of dg the vortices are pinned exactly to 
the plaquettes, while for non-integer values, as illustrated in Figure [2]|^a)-(d), they are 
interpreted being located at some intermediate position. As dg 0, the vortices are 
brought to the same plaquette. We refer to this as fusing the vortices. 

Figure [3]^b) shows the spectral evolution as a function of the vortex separation dg. 
As observed in [22l 144] , in the presence of a pair of vortices the energy of the lowest 
lying mode, S*^" = Ef% converges to zero with vortex separation. The finite energy at 
small separations is interpreted as being due to short-range vortex-vortex interactions. 
Higher energy modes are insensitive to it and remain at constant energy A/ = £"2" 
that corresponds to the spectral gap. The continuous nature of the transport, however, 
reveals the oscillations in the energy E'^" that have been qualitatively predicted in the 
context of p-wave superconductors [24]. Due to particle- hole symmetry, the oscillating 
modes come in ±£^^= pairs, but the crossings are genuine. In general, the oscillating 
energy can be written as 

£'^' = AfCos{—^)e «, (5) 
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Figure 4. (a) A plot of ln(f for J = 1 and K = 0.05. A linear fit on the peaks 
gives ^ « 2.6 (in units of ds) for the coherence length. The distance between the 
shown minima gives for the oscillation wavelength A ~ 6.0. (b) The coherence length 
^ (blue solid line) and the oscillation wavelength A (black dashed line) as functions of 
K for several J = Jx/Jz = Jyl Jz- For K < 0.1 ^ scales as K~^, while A remains a 
J dependent constant. For if > 0.1 the coherence length ^ converges to a constant, 
while A acquires vortex separation dependence (this appears as an increase due to A 
being calculated from fixed nodes at small ds for which the wavelength first increases), 
before also converging to a J independent constant. 

where non-universal constants X{J,K) > and ^{J,K) > depend on the couphngs 
and parametrize the frequency of the oscillations and the convergence of the energy, 
respectively. We are especially interested in the latter as it gives the characteristic 
coherence length for particular values of the couplings. Only when dg S> ^, one expects 
the low energy physics to be fully described by the topological Ising anyon theory. 
The parameters A and ^ can be extracted from the plot of ln(£^'^'') as shown in Figure 
Hl^a). The linear fit with negative slope confirms the exponential convergence of the 
energy with vortex separation, and distance between successive dips gives the half of 
the wavelength. By performing similar linear fits for a range of J and K, we obtain 
Figure that shows and A as functions of the physical parameters. They show very 
different behavior in the small and large K limits, which roughly occur when K < 0.1 
or K > 0.1, respectively. 

Let us consider first the K < 0.1 regime, that is experimentally the more relevant 
one. There the behavior of both the coherence length and the oscillation wavelength 
is simple. We find ^ decreasing smoothly as K~^, with the J dependence becoming 
significant only for K < 0.01. This behavior can directly be attributed to the A ~ 7^ 
scaling of the spectral gap, which implies that the coherence length scales in general as 
C, ~ A^^. For instance, when J = 1 the gap is given by A = 6y/3K (see Figure ID^b)), 
which provides an excellent approximation of the K <0.1 region in Figure Hl^b). While 
the coherence length is predominantly controlled by K, the oscillation wavelength A 
depends only on J in this regime. This follows from A being proportional to the inverse 



Interacting non-Abelian anyons as Majorana fermions in the honeycomb lattice model 9 

Fermi momentum [24]. To be precise, here it is proportional to the difference of the 
two Fermi point momenta, i.e. ^ ~ IPf ^ PfI- ■^'^^ small K the Fermi momenta are 
insensitive to the three spin term, while they move away from each other, i.e. |p^ — Pp\ 
increases, when one approaches the phase boundaries (J — )■ |) [32]. Thus the observed 
decrease in A with decreasing J is consistent with this picture. 

The behavior in the i^' > 0.1 regime is dominated by the behavior of the spectral 
gap. It converges to a J dependent constant value somewhere in the range 0.1 < K < 0.2 
and can not be increased further by increasing K [22] • When this occurs ^ becomes 
also a constant on the order of the lattice spacing. On the other hand, the oscillation 
wavelength A enters an intermediate regime where it becomes first dependent on the 
vortex separation dg. In Figure Hj^b), where the wavelength is calculated from the small 
dg behavior, this corresponds to the regime where A first increases and then decreases. 
When K is increased further, the large dg behavior takes over and A converges to the 
value of two lattice constants. This means that for sufficiently large K the oscillations 
are essentially absent when vortices are pinned to plaquettes (integer values of dg). This 
behavior is in agreement with the results for p-wave superconductors, where changes 
in the oscillations are known to occur both in the vicinity of phase transitions and for 
suffiently large gaps [211 [25]. Even though the gap saturates in the honeycomb lattice 
model, the topological phase remains robust even for K ^ J [12]. We interpret the 
changes in ^ and A being due to the specific microscopies of the honeycomb lattice model 
that can modify the non-universal properties of the topological phase. While both the 
small and large K regimes support the same non-Abelian anyons, they can therefore 
exhibit very different microscopic signatures depending on the physical parameters. 
Understanding the microscopies that control these signatures is relevant, for instance, 
to anyon-anyon interaction driven phase transitions |29j, as well as to any topological 
quantum computing schemes where read-out of information is based on detecting the 
energy shifts of the topological states [2^ . 

3.1.1. Rotational anisotropy and transversal transport To characterize the fusion mode 
energy E'^" for arbitrary relative vortex locations, we consider two further transport 
processes. First, let us consider the spectral evolution for the transport shown in Figure 

for various orientations of the lattice. If we assume the vortex at the lower left corner 
to be located at the origin, we find energy behaving identically in six sectors bounded by 
the directions ^{jo-x + n^,), -^{jo-x — 2ny) and -^{2yIx — 2ny). Figures [5l^b) and[5]^c) show 
that within these sectors one can define parallel directions along which the oscillations 
have always the same wavelength with the nodes coinciding. The only difference is the 
overall exponential damping with the vortex separation when the transport is carried 
out further away from the origin. This sector structure contrasts with the rotational 
isotropy in p-wave superconductors [15]. The difference is again in the microscopies, 
which this time derives from to the underlying honeycomb lattice. 

Finally, we consider the fusion mode energy when vortices are transported 
transversally with respect to each other as illustrated in Figure [6]|^a). As the oscillations 
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Figure 5. (a) Transportation of a vortex at separation x away from the sector 
boundary (black dashed line), (b) When a vortex is fixed at the origin located at the 
lower left corner and the other is transported along the arrows, the fusion mode energy 
behaves identically along the arrows of same colour. Parallel arrows form altogether 
six sectors with identical behavior. The picture shows parts of two sectors that are 
separated by a boundary in direction lix + ny. (c) When vortices are moved in parallel 
directions as shown in (b), the energy oscillates with the same wavelength. The plot 
colours correspond to the arrows of same colour in (a). 




Figure 6. (a) A parametrization for transversal transport. Starting with a 
configuration of two vortices separated linearly by x plaquettes, we define the 
transversal direction to be — n^; + riy. denotes the distance from the initial 
configuration, (b) For an even x the energy splitting may oscillate, but one always 
arrives at the same sign of the splitting. For an odd x the sign always changes such 
that the energy splitting is exactly zero at the sector boundary. 

are known to behave identically along the sector centers, one expects the energy to be 
reflection symmetric with respect to the boundaries. This is indeed the shown 
in Figure Mjo) for various distances x from the origin, but we also find that the sign of 
the energy may change. For an even x the energy may oscillate, but one always arrives 
at the same sign, while for an odd x the sign always changes such that the energy 
is exactly zero at the sector boundary. This implies that one can arrive at different 
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Figure 7. (a) The low-energy spectrum of the two vortex sector as a function of dg 
for J = 1 and K — 0.05. All energies are with respect to Eq, the ground state energy 
of the vortex- free sector. The dashed black line is the fermion gap Af in the absence 
of vortices. The solid (squares) and dashed (circles) red lines are the exponentially 
degenerate two vortex sectors states with energies £'q° and Eq" + , respectively, 
that correspond to an unoccupied and an occupied fusion mode. The markers denote 
positions when the vortices are located exactly at plaquettes. (b) The fermion gap A/, 
the vortex gap and the lattice potential V as functions of K. The dashed line is 
the analytic solution for the fermion gap [10] , that in the presence of vortices vanishes 
at K — only in the thermodynamic limit. For K > 0.2 the fermion gap becomes a 
constant A/ « 2.0. 



conclusions on the energy depending on the path chosen to transport the vortices. As 
all physical vortex states have to be gauge invariant, i.e. not depend on the path [TO] , 
this sets limits on simulating vortex transport through manipulating the couplings. Only 
tuning the couplings along paths without an ambiguity simulate the evolution of the 
physical states. This can be achieved by considering always the shortest path, which 
corresponds to tuning only the Jy and couplings. The restriction to two out of the 
three link types arises also naturally in the fermionization techniques that do not require 
additional gauge symmetrization [361 EH] • 

3.2. The low energy spectrum and the fusion degrees of freedom 

It has been argued in [22] that the modes having the energy behavior (j5]), to be referred 
to here as fusion modes, should be identified with the anyonic fusion degrees of freedom. 
We explicitly verify this by plotting in Figure [7|^a) some of the lowest lying states in 
the two vortex sector with respect to the lowest lying states in the vortex-free sector 
{wp = 1 on all plaquettes). At large dg there are two exponentially degenerate states 
with energies E^" and E^" + S'^" that differ by the occupation of the fusion mode. Here 
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Eq" denotes the energy of the state with unoccupied fusion mode. We define it by 



as the state with the lowest energy when vortices are nearby. It follows from this 
definition, that when the vortices are brought closer and finally fused, the energy of the 
state with the unoccupied fusion mode first oscillates and then evolves smoothly to £"0, 
the ground state of the vortex-free sectors. On the other hand, the energy of the state 
with the the occupied fusion mode becomes Aj, the first excited free fermion state on 
the vortex- free sector. In the language of topological field theories P,[in], the vortex-free 
ground state, the vortices and the fermions carry a topological quantum numbers 1, a 
and ip, respectively. They satisfy the fusion rule a x a = 1 + ip, which tells that there 
are two alternatives how a pair of vortices can behave when they are fused: they can 
either annihilate to the vacuum or leave behind a fermion. Figure [7](a) demonstrates 
that this fusion degree of freedom is exactly that of the occupation of the fusion mode^. 
We also observe that when the vortices are fused, the vacuum channel is energetically 
favoured. This contrasts with the behavior of quasiholes in the Moore-Read state, for 
which one finds the fusion channel to be energetically favoured [23]. These different 
results highlight the different microscopies of the models, that can dramatically affect 
the non-universal behavior of the topological phases. 

While the honeycomb lattice model does not accommodate analytic treatment of the 
fusion modes, their origin can be understood in the context of p-wave superconductors to 
which the honeycomb lattice model can be mapped |35| [36 t 138] . There one can explicitly 
show that vortices bind unpaired massless Majorana fermions [9], whose wave functions 
exhibit exponentially damped oscillating tails ^45j. The width of the wave functions is 
controlled by the fermion gap Aj that can be viewed as a height of a potential well 
confining them to the vortex cores 0. When two vortices are nearby each other, the 
overlap of these tails results in a finite tunneling amplitude between vortex cores, which 
in turn lifts the degeneracy and gives rise to an oscillating finite energy [24]. Our 
results of the oscillating fusion mode energies are consistent with this picture. While 
the underlying lattice makes it hard to visualize the oscillations in the wavef unctions, 
by taking appropriate linear combinations of the fusion modes one can find Majorana 
modes that have support only on the sites around an isolated vortex [16]. In Section 4 
we will employ this picture of localized Majorana modes to construct an effective lattice 
model for studying degeneracy lifting in many vortex systems. 

II Strictly speaking there is no ground state degeneracy in the two vortex sector due to fermionic parity 
conservation. However, the same study could have been carried out with another vortex pair in the 
system, in which case it is possible to define degenerate states of same fermionic parity. As long as 
these additional vortices were far away from the other, this would not affect the spectrum and the 
presented analysis would apply identically. 

% This is also in agreement with the general theory that says that the fusion channel degeneracy lifting 
can always be understood as virtual exchange of a fermionic ijj quasiparticle |21j . As Ay gives their 
minimal energy, such processes become less probable as the fermions become more massive. 




(6) 



i=2 
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3.3. Energy gaps and stability 

When the vortices are well separated, the topological low energy spectrum of Figure 
[Tl^a) can be characterized by the following energy gaps: 

A/ =^2% (7) 
A2. = lim E^,' - Eo, (8) 

dg—^oo 

V = lim E^'^' -E^'. (9) 

Af and A2,, are the fermion and vortex gaps, respectively, that describe the minimal 
energy required to excite a free fermion mode and a pair of non-interactiong vortices. V 
is a local potential that favours the vortices to be located at the plaquettes. To study 
their physical parameter dependence, we plot them in Figure ©(b) as functions of the 
effective magnetic field K. We observe a crossing in the magnitudes of the energy gaps 
such that for K < 0.1 (K > 0.1) there holds Af < A21, {Af > A2v). This means that 
while for all values of K the vacuum channel is energetically most favoured, for K < 0.1 
the energy of the system can also be lowered by fusing the vortices into fermions. For 
large effective magnetic field {K > 0.2) the energy gap becomes a constant, the exact 
value depending only on J [22\ {Af ^ 2.0 for J = 1), and can no longer be made larger 
by increasing K. 

The periodic "hopping" of the energies as a function of dg has not been observed 
before. The minima always occur for integer values of dg, i.e. when vortices are pinned to 
plaquettes, whereas the maxima occur at ds/2, i-e. when the transported vortex occupies 
a composite plaquette twice the size of a regular plaquette. We interpret this energy 
required to move a vortex to an adjacent plaquette as a local potential of magnitude V, 
defined asymptotically by (I9l). As shown in Figure [7](b), it varies only slightly with K 
and can thus be effectively treated as constant for a given uniform coupling configuration 
J. Therefore, the vortices can be viewed as living on a uniform periodic background 
potential whose minima coincide with the sites of the dual triangular lattice. 

Together the energy gaps Af and A21, and the potential V give a measure of the 
stability of the topological low energy theory against local perturbations. When the 
perturbations are weaker than Af and A2,,, the spontaneous creation of both fermion 
and vortex excitations is suppressed. When they are also weaker than V, any existing 
vortices in the system will remain stationary. A / describes the stability of the topological 
phase itself, whereas A2^ and V describe stability of the vortex sectors within it. 

4. Degeneracy lifting in many vortex systems 

In this section we study degeneracy lifting in many vortex systems. Employing an 
effective lattice model of Majorana fermions, we show that the hybridized spectrum 
of the fusion modes can be understood in terms of bi-partite interactions between the 
vortices. 
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Figure 8. (a) A four vortex configuration withi one vortex separated by ds from the 
other three, (b) One of the two fusion modes oscillates and converges quickly to zero 
energy, while the other remains localized at finite energy on the group of three vortices. 

4.I. Free Majorana fermions on a lattice 

In the previous section we analyzed the spectral evolution in two vortex systems. A 
simple generalization is to consider what happens when a single vortex is dragged away 
from a bunch of vortices. When this is performed for the four vortex configuration 
shown in Figure [Hl^a), one still finds qualitatively similar behavior. Figure El^b) shows 
that one of the fusion modes to oscillates and converges to zero energy, while the other 
remains at finite energy. In other words, one of the two fusion modes is delocalized while 
the other remains localized at the three vortices. However, comparing to Figure Et^b), 
the energy of both fusion modes, including the amplitude and the wavelength of the 
oscillations, is modified due to the presence of additional vortices. As more vortices are 
added, the fusion modes will in general hybridize a new energy band within the spectral 

gap m- 

The degeneracy lifting due to tunneling Majorana fermions suggests that this energy 
band can be modelled by Majorana fermions tunneling on a lattice defined by the vortex 
locations [31]. Such a model possesses two properties that we have already observed 
on the honeycomb lattice model: it has particle-hole symmetry giving rise to double 
spectrum and it exhibits exact zero energy states for an odd number of sites. The most 
general Hamiltonian one can write is 

H = tJ2Yl ^'47.7,' (10) 
I \i-j\=i 

where ti are tunneling amplitudes between sites of separation /, s-^ are local Z2 gauge 
degrees of freedom on link (ij) and 7^ are Majorana fermions satisfying 7^! = 7j and 
{ji, 7j} = 26ij. The Hamiltonian can be defined on arbitrary lattice geometries, which in 
our case is fixed by the positioning of the vortices. As the interactions in the honeycomb 
lattice are exponentially suppressed, we simplify the effective model by considering only 
the three shortest range interactions of range / = 1, a/3, 2. We refer to them as ti, t^ 
and t2 interactions, respectively, that couple sites as illustrated in Figures |9]^c) and|9]^f). 
The configurations = {s-j } are called gauge degrees of freedom, because the spectrum 




Interacting non-Abelian anyons as Majorana fermions in the honeycomb lattice model 15 



does not depend on them explicitly. Instead, the model depends on the introduced 
effective flux, which we define on plaquette [ijk), i, j and k being the three corners of 
a triangular plaquette enumerated counter-clockwise, by 

$ijfc = -iln(iSijSjfcSfci). (11) 

This evaluates always to either ±|. The set of fiuxes on all plaquettes is referred to as 
a flux sector. 

The free parameters of the model are the couplings ti and the fiux sectors. Our 
objective is to study how they are to be fixed such that the spectrum of the Majorana 
model will approximate the spectrum of the fusion modes. To compare quantitatively 
the spectrum Sn = {~En/2, —En/2-1, ■ ■ ■ , En/2-1, En/2) of the fusion modes from an n- 
vortex system to the spectrum e„ = (—en/2, —^n/2-1, ■ ■ ■ , ^n/2-1, en/2) of a corresponding 
n-site effective model, we introduce the deviation measure 

F[£n,en]=-\£„-en\. (12) 

For F = the spectra match exactly. To systematically study the many vortex 
systems, we restrict to considering homogenous vortex configurations, denoted by 
^N^xNyi consisting of vortices in direction n^. and A^^^ vortices in direction n^, all 
separated by at least d links. In particular, we will consider chain (A^j, = 1) and ladder 
{Ny = 2) configurations that are illustrated in Figures |9]^a) and Mjo), and Figures 
Ht^d) and M^e), respectively. Configurations that differ only by the sparsity d have 
topologically the same effective description. For instance, all four vortex chains are 
described by an effective Majorana model living on the four site lattice shown in Figure 

Etc). 

Our strategy for studying the many vortex degeneracy lifting is as follows. We set 
again globally Jx = Jy = Jz = ^, but as before leave the effective magnetic field strength 
as a free parameter. Considering then vortex chain and ladder configurations, we 
find the fusion mode spectra S{K) for the range < < 0.3. Since the fusion mode 
energy S'^" in two vortex systems can be understood as arising due to tunneling between 
vortex cores, its magnitude at a fixed vortex separation dg = I should correlate with the 
tunneling amplitudes between sites of separation /. We assume the simplest possible 
correlation and use the ansatz 

ti{K)=£'^{K), (13) 

for the effective model couplings corresponding to Cfj^^^^ vortex configuration. For 
the fiux sectors, on the other hand, there is no obvious a priori correlation with the 
underlying model. We assume only that topologically equivalent effective models, i.e. 
ones differing only by the sparsity (i, should be defined on the same fiux sector. To 
find the one providing the correct effective description, we need to consider the effective 
spectrum e{K) = e[ti{K)] over all the possible different fiux sectors. 
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(d) (e) (f) 

Figure 9. The (a) dense C^^i and (b) sparse C|xi chain vortex configurations 
have both (c) the topologically same effective Majorana model description. Likewise, 
(d) and (e) are six vortex ladder configurations C|x2 and C|x2, respectively, both 
described by the same effective model shown in (f). The solid, dashed and dotted lines 
denote ti, t2 and interactions, respectively. To evaluate the flux per plaquette, we 
use counter-clockwise convention to assign the overall i (— i) factor per link when a 
Majorana fermion hops along (against) the shown orientation. 

4-1 .1. Vortex chains Let us consider first vortex chains where the t^ interactions are 
naturally absent due to the vortex arrangement. Figures [TUT a) and lTUT b) show the fusion 
mode and the effective Majorana model spectra with nearest neighbour interactions only 
for the dense Cl^-^ and sparse Cf^i four vortex chains, respectively. In the absence of 
also ^2 interactions, there are no non-trivial plaquettes and thus the spectrum does not 
depend on the flux sector. We find that the ti interactions alone provide an excellent 
approximation, which, as quantitatively shown in Figures [TUT c) and fTUT d). corresponds 
to at least F < 0.02 and gets better with increasing K. This is in agreement with longer 
range tunneling being more suppressed for larger K. We note also that in Figure [TOlfb) 
the degeneracy point in the fusion mode spectrum S{K) and the zero energy point of 
the vortex- vortex interaction induced fusion mode energy coincide exactly at K ~ 0.12. 
This strongly suggests that the nearest neighbour tunneling is primarily responsible for 
the degeneracy lifting in vortex chains. 

When the ^2 interactions are switched on. Figures [TUlfc) and [TUlfd) show that the 
approximation can be improved for i^' < 0.1 only if one chooses the fluxes alternating 
between | and — |. For instance, for the C^^i configuration shown in Figure [9]^c) one 
should choose [$123, '^'243] = ilf;"!]- the other hand, for K > 0.1 we find that 
the inclusion of the ^2 interactions, regardless of the flux sector, has either neglible 
or detrimental effect on the approximation. The distinct behavior at the different K 
regimes can again be traced back to the K controlled localisation of the Majorana modes. 
For the experimentally interesting small K regime the longer range interactions are more 
relevant and should be included in the effective model. On the other hand, for large K 
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Figure 10. The fusion mode and effective model spectra for (a) the dense Cj^i 
and (b) the sparse Cf^i chain configurations. The red dashed hnes denote fusion 
mode spectra £, the black dash dotted lines the bi-partite fusion mode energies S"^' 
corresponding to ti interactions, and the squares and circles the effective model spectra 
with (on the alternating flux sector) and without the t2 interactions, respectively. The 
deviations F for (c) the dense chains Cj^^j^ and (d) the sparse chains Cff^^ for different 
flux sectors. The circles correspond to ti effective model only, while the crosses and 
squares correspond to uniform and alternating flux configurations, respectively. The 
approximation by the effective model compared to the pure ti model can only be 
improved when employing the alternating flux sector. The region K < 0.03 is omitted 
due to finite size effects. 

they are either neghgible due to ^2 ^ (sparse d = 2 chains) or the bi-partite sphtting 
f lT3|) overestimates their strength due to collective effects (tight d = 1 chains). As these 
observations apply also to longer vortex chains, we conclude that the ^2 interactions are 
physically relevant in the small K regime of the vortex chain systems. The necessary 
flux sector is the one having ±| flux alternating on the plaquettes consisting of two 
ti-links and a single t2-link. 

4- 1-2. Vortex ladders The main difference between the vortex chains and the ladders 
is the additional presence of the t^ interactions. However, let us start again with 
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Figure 11. The fusion mode and effective model spectra for (a) tfie dense 6*2x2 
and (b) the sparse C|x2 vortex ladder configurations. The red dashed lines denote 
the fusion mode spectra £, the black dash dotted lines the fusion mode energies E"^' 
corresponding to ti interactions, and the circles and squares the ti effective model 
spectra on the uniform and the alternating flux sectors, respectively. The deviations 
F for (c) the dense C^x2 ^i^*^ (d) the sparse C^x2 ladders for different flux sectors. 
The circles and squares denote flux sectors as above. The region K < 0.03 is omitted 
due to finite size effects. 



the nearest neighbour effective model which has now two distinct ffux sectors: the 
uniform sector with $125 = $256 = ^235 = '^'345 or the alternating sector with 
$126 = $235 = -$256 = -$345 (see Figurel^f )). Figures [111(a) and [111(b) for the smallest 
dense and sparse ladders, respectively, show that a qualitatively good fit can be obtained 
in the uniform flux sector, whereas in the non-uniform sector the spectrum is gapless. 
This is reflected in the deviation measures, shown in Figures [TlTc) and [T]T d). where the 
uniform flux sector clearly provides a better approximation than the non-uniform one. 
Therefore, we conclude that having a |-flux on all plaquettes consisting of ti links only 
is another necessary condition for the effective Majorana model to approximate the full 
model. 

In general the approximation by the nearest neighbour model is not as good for 
the vortex ladders as it was for the vortex chains. However, this error is in general 



Interacting non-Abelian anyons as Majorana fermions in the honeycomb lattice model 19 




(a) (b) (c) 

Figure 12. Flux on the effective model arises from the tt-Aux vortices of the 
underlying model. When the flux of a vortex is assumed to be uniformly spread across 
the plaquette, the flux through the shown effective model plaquettes is given by (a) 

6 6 6 2 ' ^ u vl; ^ 12 ^ 12 3 2 ' 



systematic and can be compensated by scaling the tunneling amplitudes as ti = aS 
for some K indepedent constant a. While the systematic study of such dressing 
effects is beyond the present work, one can imagine them occurring due to quantum 
interference of many Majorana wave functions. Such effects are also likely to play a 
role when the and ^2 interactions are switched on. In this case the flux needs to 
be defined also on plaquettes consisting of two ti-links and a single t^-link, and on 
plaquettes with one link of each type. For instance, in Figure Mf) these correspond to 
the fluxes [$155, $125, $245, '^'234, $356, $362] and [$136, $135, $462, $463]- We find that the 
approximation can in general be improved for some K ranges by introducing the longer 
range interactions, but unlike in the vortex chains we find no single flux sector that would 
systematically improve it in the physically relevant small K regime. We attribute this 
to the studied vortex arrays being too small for such genuinely two dimensional effects 
to become transparent. We leave their systematic study for future work |17] . 

4-2. Effective flux sectors and the underlying n-flux vortices 

We discovered that having a |-flux on every plaquette consisting of only ti-links is a 
necessary condition for the effective Majorana model to approximate the full one. Also, 
in the presence of t2 interactions, the approximation for vortex chains could be improved 
only when one imposed an alternating ±|-fiux pattern on the plaquettes consisting of 
two ti-links and a single t2-link. We postulate that these conditions are connected to the 
underlying honeycomb lattice model as follows. The anyonic vortices there are vr-fiux 
vortices. Let us assume the flux to be uniformly spread across the hexagonal plaquette 
a vortex occupies and the Majorana mode to be bound exactly at its center. Then, as 
illustrated in Figure [T2t^a), a triangular dual lattice plaquette consisting of three ti-links 
will enclose exactly ^-fiux arising collectively from the three vortices. On the other 
hand, the plaquettes consisting of a single t2 link and two ti links, as illustrated in 
Figure [TW b). span no area and should thus have no flux on them. While trivial flux 
can not directly imposed on the Majorana model, we interpret the alternating ±|-flux 
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sector effectively providing it. 

We postulate that when writing down a general effective Majorana model for 
interacting 7r-ffux vortices, the flux on every possible effective model plaquette should be 
chosen to be consistent with the enclosed vortex ffux from the underlying microscopic 
model. For instance, as illustrated in Figure [T^c), the plaquettes consisting of a single 
t^-link and two ti-links should have also |-flux. While the concept of enclosed vortex 
flux is equivalent of the net phase in vortex arrays [31], we emphasize the conceptual 
simplicity of our prescription in relating the microscopic model and the necessary 
conditions it imposes on the effective model. The verification of this picture and the 
development of a thermodynamic Majorana model is a subject of future research [17] . 

5. Conclusions 

We have studied the interacting non-Abelian Ising anyons emerging from Kitaev's 
honeycomb lattice model. By simulating continuous vortex transport [30] , we uncovered 
the characteristic oscillations in the bi-partite degeneracy lifting and characterized their 
physical parameter dependence. While the oscillations had been discovered earlier for 
both p-wave superconductors [25] and the Moore- Read state [23], these calculations 
involve mean fields and trial wave functions, respectively. Our results demonstrate the 
oscillations in an actual microscopic model. Due to algebraic similarity, we expect them 
to be present also in the variations of the honeycomb lattice model [HI [121 lEl El US] ■ 
Furthermore, we extended the results of [22] by demonstrating unambiguously the fusion 
degrees of freedom and by obtaining the energy gaps characterizing the stability of the 
topological low energy theory. 

In the second part we studied degeneracy lifting in interacting many vortex systems. 
Employing the picture of localized Majorana modes, we wrote down a Majorana fermion 
model on a finite lattice whose sites coincide with the vortex locations. By comparing its 
spectrum to that of the fusion modes of a corresponding vortex configuration, we were 
able to find necessary conditions for the spectra to match. This amounted to finding 
relevant flux sectors of the effective model, which we interpreted as arising from the tt- 
flux vortices of the underlying full model. These results are consistent with and extend 
the earlier mean-field studies in the context of p-wave superconductors [31]. Further, 
we found that the energy splitting due to vortex-vortex interactions could be employed 
to high accuracy as the effective tunneling amplitude of the Majorana fermions. This 
confirms that bi-partite tunneling is predominantly responsible for the degeneracy lifting 
also in many vortex systems. If many body effects are present, they can be included 
as dressed or longer range couplings. Our results pave the way for constructing an 
effective model for interacting anyonic vr-flux vortices at the thermodynamic limit [47j . 
Such model can be employed to understand phase transitions due to vortex lattices, 
[32] . and test the general theory of anyon-anyon interaction driven phase transitions in 
the context of the microscopic honeycomb lattice model [29l H?] . 
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